% CSP function
% X1 is the EEG signal from target 1 while X2 is the other category
% The data format must be channel in rows and time in columns;

% m is the number of characteristics you want to retain in W;

% Wr is the spatial fiter, its inverse inv(Wr) is the common spatial
% pattern;
function Wr = common_spatial_pattern(X1,X2,m)
cov1 = zeros(size(X1,1),size(X1,1),size(X1,3));
cov2 = zeros(size(X2,1),size(X2,1),size(X2,3));
for i = 1:size(X1,3)
    tmp = X1(:,:,i);
    cov1(:,:,i) = tmp*tmp'/trace(tmp*tmp');
end
for j = 1:size(X2,3)
    tmp = X2(:,:,j);
    cov2(:,:,j) = tmp*tmp'/trace(tmp*tmp');
end
cov1 = sum(cov1,3);
cov2 = sum(cov2,3);
covTotal = cov1 + cov2;
[Ut, Lt] = eig(covTotal);
Lt = diag(Lt);
[Lt,LtIndex] = sort(Lt,'descend');
Ut = Ut(:,LtIndex);
P = diag(sqrt(1./Lt))*Ut';
S1 = P*cov1*P';
[U1, L1] = eig(S1);
L1 = diag(L1);
[~,L1Index] = sort(L1,'descend');
U1 = U1(:,L1Index);
W = U1'*P;
Wr = W([1:m (end-m+1):end],:);
end